#include "TCanvas.h"
#include "THStack.h"
#include "TList.h"
#include "TLatex.h"
#include "BasicTypes.hxx"
#include "DataSample.hxx"
#include "DrawingTools.hxx"

#include "NeutrinoUtils.h"
#include "load.h"

Bool_t fkLogy = kFALSE;


void saveHist(const TString opt, TList *lout, TCanvas *c1, Int_t &id, const TString tag, const Int_t accuml, TString tit, const Bool_t k2natX=kFALSE, const Bool_t k2natY = kFALSE)
{
  //c1->GetListOfPrimitives()->ls();
  //printf("test size %d\n",c1->GetListOfPrimitives()->GetEntries());
  c1->SetTopMargin(0.1);
  const TString tt = tit(0,tit.First(";"));
  TLatex *lt = new TLatex(0.1,0.92,tt);
  lt->SetNDC();

  TLegend * lg = (TLegend*)(c1->GetListOfPrimitives()->At(c1->GetListOfPrimitives()->GetEntries()-1));
  lg->SetName(Form("lgID%02d_%s_%d", id, tag.Data(), accuml));
  tit=tit(tit.First(";"), tit.Length());
  if(!lg){
    printf("lg null!!\n"); 
    c1->GetListOfPrimitives()->ls();
    exit(1);
  }
  lout->Add(lg);

  THStack * hsin = (THStack*)(c1->GetListOfPrimitives()->At(1));
  TList *ltmp = hsin->GetHists();
  //ltmp->ls();

  THStack * hsout = new THStack(Form("ID%02d_%s_%d", id, tag.Data(), accuml),tit);
  lout->Add(hsout);

  for(Int_t ii=0; ii<ltmp->GetEntries(); ii++){
    TH1 * obj = (TH1*)(ltmp->At(ii));
    if(!obj){
      printf("%d null\n",ii);exit(1);
    }

    obj->SetTitle(tit);
    if(k2natX){
      NeutrinoUtils::ToNaturalScale(obj->GetXaxis());
    }
    if(k2natY){
      NeutrinoUtils::ToNaturalScale(obj->GetYaxis());
    }
    hsout->Add(obj);
    printf("%d color %d %d\n", ii, obj->GetFillColor(), obj->GetLineColor());
  }

  if(opt.Contains("box")){
    //opt+=" nostack";
  }
  hsout->Draw(opt);
  printf("%s opt %s\n", hsout->GetName(), opt.Data());

  if(k2natX){
    gPad->SetLogx();
  }
  else{
    gPad->SetLogx(0);
  }

  if(k2natY){
    gPad->SetLogy();
  }
  else{
    gPad->SetLogy(0);
  }
  
  if(fkLogy){
    gPad->SetLogy();
  }

  lg->SetX1NDC(0.85);
  lg->SetX2NDC(0.99);
  lg->SetY1NDC(0.65);
  lg->SetY2NDC(0.99);
  lg->Draw();

  lt->Draw();

  c1->Print(Form("outplot/ID%02d_%s_accumlevel%d.png", id, tag.Data(), accuml));
  
  id++;

  gPad->SetLogx(0);
  gPad->SetLogy(0);

  delete lt;
}

void plothighland2(const string fin, const Int_t accuml = 5)
{
  load();
  
  DrawingTools draw(fin);
  DataSample dsm(fin);
 
  TTree * tt = (TTree*) gDirectory->Get("default");
 
  TCanvas *c1=new TCanvas("c1","",600,400); 

  const TString mpt("selmu_mom*sin(selmu_theta)/1e3");
  tt->SetAlias("truemu_trueE", "sqrt(truemu_truemom*truemu_truemom+105.65837*105.65837)");
  const TString mQ2("( 2*( nu_trueE*truemu_trueE - truemu_truemom*nu_trueE*truemu_truecostheta ) - 105.65837*105.65837 )/1e6");
  const TString mTruePt("truemu_truemom*sqrt(1-truemu_truecostheta*truemu_truecostheta)/1e3");
  //===============

  //tt->SetAlias("mu_mass", "105.65837");  //fail
  tt->SetAlias("linearQ2",mQ2);
  tt->SetAlias("selmu_pt", mpt);
  tt->SetAlias("truemu_truept", mTruePt);

  tt->SetAlias("lpt",Form("TMath::Log10(%s)", mpt.Data()));
  tt->SetAlias("lTruePt",Form("TMath::Log10(fabs(%s))", mTruePt.Data()));
  tt->SetAlias("lQ2", Form("TMath::Log10(%s)", mQ2.Data()));
  tt->SetAlias("lmom","TMath::Log10(selmu_mom/1e3)");  

  string cut0(Form("accum_level > %d", accuml));
  string lev0("particle");


  ///---
  const Int_t enbin = 60;
  const Double_t emin = -0.1;
  const Double_t emax = 5;

  const Double_t lenbin = 25;
  const Double_t lemin = TMath::Log10(0.1);
  const Double_t lemax = TMath::Log10(20);
 
  const Double_t ldnbin = 25;
  const Double_t ldmin = -1;
  const Double_t ldmax = 1;
  //---

  TList * lout = new TList;
  //lout->SetOwner(kTRUE);
  Int_t hid=0;
  TString tag;

  fkLogy = kTRUE;
  draw.Draw(dsm,"nu_trueE/1e3", enbin, emin, emax, lev0, cut0); tag="nuE"; saveHist("hist",lout, c1, hid, tag, accuml,";E_{#nu, sim} (GeV);Counts"); 
  draw.Draw(dsm,"selmu_mom/1e3", enbin, emin, emax, lev0, cut0); tag="mom"; saveHist("hist",lout, c1, hid, tag, accuml,";p_{sel.#mu, rec} (GeV/c);Counts");
  draw.Draw(dsm,"selmu_pt", enbin, emin, emax, lev0, cut0); tag="pt"; saveHist("hist",lout, c1, hid, tag, accuml,";p_{t}^{sel.#mu, rec} (GeV/c);Counts");
  draw.Draw(dsm, mTruePt.Data(), lenbin, lemin, lemax, lev0, cut0); tag="truept"; saveHist("hist",lout, c1, hid, tag, accuml,";p_{t}^{true #mu, sim} (GeV/c);Counts");
  draw.Draw(dsm,"selmu_theta*TMath::RadToDeg()", 45, 0,90, lev0, cut0); tag="theta"; saveHist("hist",lout, c1, hid, tag, accuml,";#theta_{sel.#mu, rec} (deg.);Counts"); 

  fkLogy = kFALSE;
  draw.Draw(dsm, "selmu_mom/truemu_truemom-1:lQ2", lenbin, lemin, lemax, ldnbin, ldmin, ldmax, lev0, cut0,"box"); tag="momRes"; saveHist("box",lout, c1, hid, tag, accuml,";Q^{2}_{sim} (GeV^{2}/c^{4});p_{sel.#mu, rec}/p_{true #mu, sim}-1",kTRUE);
  draw.Draw(dsm, "cos(selmu_theta)/truemu_truecostheta-1:lQ2", lenbin, lemin, lemax, ldnbin, ldmin, ldmax, lev0, cut0,"box"); tag="costhetaRes"; saveHist("box",lout, c1, hid, tag, accuml, ";Q^{2}_{sim} (GeV^{2}/c^{4});cos#theta_{sel.#mu, rec}/cos#theta_{true #mu, sim}-1", kTRUE); 
  draw.Draw(dsm, "sin(selmu_theta)/sqrt(1-truemu_truecostheta*truemu_truecostheta)-1:lQ2", lenbin, lemin, lemax, ldnbin, ldmin, ldmax, lev0, cut0,"box"); tag="sinthetaRes"; saveHist("box",lout, c1, hid, tag, accuml, ";Q^{2}_{sim} (GeV^{2}/c^{4});sin#theta_{sel.#mu, rec}/sin#theta_{true #mu, sim}-1", kTRUE); 
  draw.Draw(dsm, "lpt/lTruePt-1:lQ2", lenbin, lemin, lemax, ldnbin, ldmin, ldmax, lev0, cut0,"box"); tag="ptRes"; saveHist("box",lout, c1, hid, tag, accuml, ";Q^{2}_{sim} (GeV^{2}/c^{4});p_{t}^{sel.#mu, rec}/p_{t}^{true #mu, sim}-1",kTRUE); 
  draw.Draw(dsm, "selmu_pt:truemu_truept", enbin, emin, emax, enbin, emin, emax, lev0, cut0,"box"); tag="ptCor"; saveHist("box",lout, c1, hid, tag, accuml, ";p_{t}^{true #mu, sim} (GeV/c);p_{t}^{sel.#mu, rec} (GeV/c)");
  draw.Draw(dsm, "lpt:lTruePt", lenbin, lemin, lemax, lenbin, lemin, lemax, lev0, cut0,"box"); tag="llptCor"; saveHist("box",lout, c1, hid, tag, accuml, ";p_{t}^{true #mu, sim} (GeV/c);p_{t}^{sel.#mu, rec} (GeV/c)",kTRUE, kTRUE); 

  fkLogy = kTRUE;
  draw.Draw(dsm,"lQ2", lenbin, lemin, lemax, lev0, cut0+"&& linearQ2>0");  tag="lq2"; saveHist("hist",lout, c1, hid, tag, accuml,";Q^{2}_{sim} (GeV^{2}/c^{4});Counts",kTRUE); 

  fkLogy = kFALSE;
  draw.Draw(dsm,"lpt:lQ2", lenbin, lemin, lemax, lenbin, lemin, lemax, lev0, cut0+"&& linearQ2>0", "box"); tag="lptlq2"; saveHist("box",lout, c1, hid, tag, accuml,";Q^{2}_{sim} (GeV^{2}/c^{4});p_{t}^{sel.#mu, rec} (GeV/c)",kTRUE,kTRUE); 
  draw.Draw(dsm,"lmom:lQ2", lenbin, lemin, lemax, lenbin, lemin, lemax, lev0, cut0+"&& linearQ2>0", "box"); tag="lmomlq2"; saveHist("box",lout, c1, hid, tag, accuml,";Q^{2}_{sim} (GeV^{2}/c^{4});p_{sel.#mu, rec} (GeV/c)",kTRUE,kTRUE); 

  //---

  string cut1 = "&& linearQ2 <1";
  string cut2 = "&& linearQ2 > 1 && linearQ2 <2";
  string cut3 = "&& linearQ2 > 2";

  fkLogy = kTRUE;

  draw.Draw(dsm,"lpt", lenbin, lemin, lemax, lev0, cut0+cut1); tag="lptcut1"; saveHist("hist",lout, c1, hid, tag, accuml,cut1+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 
  draw.Draw(dsm,"lpt", lenbin, lemin, lemax, lev0, cut0+cut2); tag="lptcut2"; saveHist("hist",lout, c1, hid, tag, accuml,cut2+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 
  draw.Draw(dsm,"lpt", lenbin, lemin, lemax, lev0, cut0+cut3); tag="lptcut3"; saveHist("hist",lout, c1, hid, tag, accuml,cut3+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 

  draw.Draw(dsm,"lmom", lenbin, lemin, lemax, lev0, cut0+cut1); tag="lmomcut1"; saveHist("hist",lout, c1, hid, tag, accuml,cut1+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 
  draw.Draw(dsm,"lmom", lenbin, lemin, lemax, lev0, cut0+cut2); tag="lmomcut2"; saveHist("hist",lout, c1, hid, tag, accuml,cut2+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 
  draw.Draw(dsm,"lmom", lenbin, lemin, lemax, lev0, cut0+cut3); tag="lmomcut3"; saveHist("hist",lout, c1, hid, tag, accuml,cut3+";p_{t}^{sel.#mu, rec} (GeV/c);Counts",kTRUE); 

  TFile *fout=new TFile("OutList.root","recreate");
  lout->Write();
  fout->Save();
  fout->Close();
  delete fout;
}


//rm -f outplot/*.png; r -q -b -l 'plothighland2.C ("test.root")'; rm -f transport/*; cp outplot/ID*.png transport/
